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Abstract 

In  this  paper  we  describe  a new  phase-field  model  to  describe  isothermal  phase  tran- 
sitions between  ideal  binary  alloy  solution  phases.  Governing  equations  are  developed 
for  the  temporal  and  spatial  variation  of  the  phase  field,  which  identifies  the  local  state 
or  phase,  and  for  the  composition.  An  asymptotic  analysis  as  the  gradient  energy  co- 
efficient of  the  phase-field  becomes  small  shows  that  our  model  recovers  classical  sharp 
interface  models  of  this  situation  when  the  interfacial  layers  are  thin,  and  we  relate  the 
parameters  appearing  in  the  phase-field  model  to  material  and  growth  parameters  in 
real  systems.  We  identify  three  stages  of  temporal  evolution  for  the  governing  equa- 
tions; the  first  corresponds  to  interfacial  genesis  which  occurs  very  rapidly;  the  second 
to  interfacial  motion  controlled  by  diffusion  and  the  local  energy  difference  across  the 
interface;  the  last  takes  place  on  a long  time  scale  in  which  curvature  effects  are  im- 
portant, and  corresponds  to  Ostwald  ripening.  We  also  present  results  of  numerical 
calculations. 


Permcinent  address:  School  Of  Mathematics,  University  Walk,  Bristol,  BS8  ITW,  UK. 


1 Introduction 


Classical  macroscopic  models  of  phase  transitions  describe  the  interface  between  regions  of 
different  phase  as  a surface  of  zero  thickness.  The  governing  equations  for  thermodynamic 
variables,  such  as  temperature  and  composition,  are  formulated  in  each  phase  independently, 
based  on  conservation  principles  and  quantitatively  verified  phenomenological  laws  that  re- 
late fluxes  to  gradients.  The  boundary  conditions  at  the  interface  are  of  two  types:  those 
which  represent  conservation  laws  when  there  are  discontinuities  in  the  thermodynamic  vari- 
ables or  their  gradients,  and  those  which  deal  with  the  values  of  the  thermodynamic  variables 
at  the  interface  and  are  based  on  modeling  of  the  interface  on  an  atomic  scale.  This  approach 
gives  rise  to  the  formulation  of  a free  boundary  problem  which  provides  a difficult  mathe- 
matical setting  in  which  only  phase  changes  with  simple  configurations  have  been  rigorously 
analysed  mathematically.  Because  these  models  have  been  used  for  many  years,  it  is  clear 
from  the  outset  what  physical  mechanisms  are  incorporated  into  them,  and  comparison  with 
carefully  controlled  experiments  have  been  performed,  e.g.  the  dendritic  growth  experiments 
of  Huang  and  GHcksman  [1]  and  the  grain  boundary  groove  measurements  of  Hardy  [2]. 

An  alternative  technique  for  investigating  systems  involving  a phase  transition  involves 
the  construction  of  a Cahn- Hilliard  or  Landau- Ginzberg  free  energy  functional,  which  unlike 
the  classical  model  treats  the  system  as  a whole.  This  approach  has  its  roots  in  statistical 
physics  (see  e.g.  Landau  and  Khalatinikov  [3]).  For  this  purpose,  a phase  field,  (^(x,  t),  is 
postulated  to  exist,  which  characterises  the  phase  of  the  system  at  each  point  in  space  and 
time.  It  is  assumed  that  the  Helmholtz  free  energy  ...],  is  a functional  of  the  phase 
field,  as  well  as  any  other  thermodynamic  variables  (such  as  temperature  and  composition 
which  are  denoted  here  by  ellipsis),  and  their  gradients,  in  the  following  way: 


•••]  = _(  [M,  ■ • ■)  + + . . .]  dn, 


(1) 


where  Q is  the  region  occupied  by  the  system,  and  /(<^, ...)  is  the  Helmholtz  free  energy 
density  for  a phase  with  no  gradients.  Its  dependence  on  usually  has  a “double- well”  form. 
Here  we  assume  that  the  only  contribution  to  the  energy  functional  from  local  gradients  is 
that  of  the  phase  field;  a more  general  model  might  also  include  contributions  due  to  the 
gradients  of  the  other  thermodynamic  variables.  The  phase  field  is  then  assumed  to  evolve 
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as: 

^oc£(^),  (2) 

where  £ is  some  partial  differential  operator,  such  that  £(0)  = 0.  This  equation  is  then 
supplemented  by  partial  differential  equations  for  the  other  thermodynamic  variables. 

In  some  situations  the  composition  naturally  plays  the  role  of  the  phase  field,  in  which  case 
the  double-well  in  the  free  energy  occurs  with  respect  to  composition.  Cahn  and  HiUiard 
[4,  5,  6]  have  used  this  approach  to  model  interfacial  energies,  nucleation  and  spinodal 
decomposition  in  a binary  alloy.  Also  Langer  and  Sekerka  [7]  have  modeled  the  motion  of  a 
planar  interface  using  this  approach. 

More  generally,  various  models  that  employ  these  ideas  are  reviewed  by  Halperin,  Ho- 
henburg  and  Ma  [8],  particularly  in  regard  to  the  study  of  critical  phenomena.  The  Model 
C given  by  Halperin  et  al.  has  been  adapted  by  Langer  [9],  Fix  [10],  and  most  prolifically 
by  Caginalp  [11],  to  derive  the  so-called  “phase-field  model”  of  solidification  which  describes 
the  phase  change  of  a pure  material.  Caginalp  has  studied  this  model,  and  its  variations 

[12,  13],  extensively.  In  this  model  the  phase  field  is  required  to  evolve  according  to: 

= = <!>")  + (3) 

where  a,  a,  and  ( are  non-dimensional  constants.  A non-dimensional  temperature,  T,  satis- 
fies the  heat  equation,  modified  to  allow  for  latent  heat  production  as  a source  term  propor- 
tional to  <^t: 

Tt  + 1^.  = KV^T,  (4) 

where  I is  the  latent  heat  production.  The  Helmholtz  free  energy  density  used  by  Caginalp 
is: 

= (5) 

where  T = 0 is  the  phase  transition  temperature. 

It  has  emerged  from  study  of  this  model  that  quaHtatively  it  exhibits  features  common  to 
soHdification  of  a pure  material.  Numerical  calculations  by  Smith  [14]  and  by  Socolowski  and 
Caginalp  [15],  and  for  a similar  model  by  Kobayashi  [16],  show  such  features  as  breakdown  of 
planar  and  circular  interfaces  to  cellular  structures,  as  well  as  the  formation  of  dendrite-like 
structures,  inclusion  of  liquid  pockets,  and  Ostwald  ripening  or  coarsening  behaviour. 
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Caginalp  [12]  has  shown  in  various  distinguished  limits,  in  which  ^ ^ 0,  that  various 
forms  of  the  classical  Stefan  problem  may  be  recovered,  in  which  the  interface  is  taken  to  be 
“sharp”  i.e.  modeled  by  a surface.  In  this  limit,  there  are  thin  layers  within  Q,  of  thickness 
0{i)  in  which  the  phase  field  rapidly  changes.  These  are  interpreted  as  representing  diffuse 
interfaces.  From  this  analysis  it  transpires  that  in  some  limits,  the  interfacial  dynamics 
involve  curvature  effects  corresponding  to  the  Gibbs-Thomson  interfacial  surface  energy  effect 
as  well  as  kinetic  effects.  In  other  hmits,  it  is  also  possible  to  recover  the  classical  Hele-Shaw 
problem.  It  is  clear  that  this  approach  can  embody  a considerable  variety  of  reafistic  physical 
effects  in  a coherent  way. 

When  applying  the  model  above  to  a definite  physical  situation  it  is  not  clear  how  to 
choose  the  values  of  parameters  in  the  phase-field  model  so  that  it  describes  the  solidification 
of  a pure  material  with  known  materials  and  growth  parameters.  Such  a choice  of  parameters 
is  essential  to  permit  comparison  with  experiment.  Also  the  free  energy  functional  is  only 
employed  in  the  formulation  of  the  kinetic  equation  for  the  phase  field  and  not  for  the 
temperature.  An  alternative  approach,  as  suggested  by  Penrose  and  Fife  [17],  is  to  construct 
an  entropy  functional  of  the  system,  and  postulate  kinetic  equations  for  (j)  and  T which 
ensure  that  the  entropy  increases  monotonically  in  time,  as  required  by  the  Second  Law  of 
Thermo  dy  nami  cs . 

The  appeal  of  phase-field  models  in  describing  phase  transitions  is  twofold: 

• It  provides  a simple,  elegant  description,  that  appears  to  embody  a rich  variety  of 
realistic  physical  growth  phenomena. 

• From  a computational  point  of  view,  it  is  relatively  simple  to  obtain  solutions.  This  is 
because  it  is  not  necessary  to  distinguish  between  the  different  phases.  Computations 
that  employ  the  classical  sharp  interface  formulation  require  the  free  boundary  to  be 
tracked  numerically,  and  the  regions  occupied  by  each  phase  are  therefore  determined 
and  dealt  with  individually.  Added  difficulties  with  the  classical  method  are  encoun- 
tered when  the  connectedness  of  the  solid  or  liquid  regions  change  in  time.  These 
considerations  often  result  in  very  difficult  and  untidy  numerical  algorithms. 

In  this  paper  we  derive  and  analyse  a new  phase-field  model  for  isothermal  phase  tran- 
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sitions  of  a solution  of  two  components.  In  an  isothermal  system  the  driving  force  for  phase 
transformation  is  developed  through  the  growth  of  a product  phase  into  a supersaturated 
parent  phase.  To  our  knowledge,  there  are  to  date  no  phase-field  models  that  deal  with  im- 
pure materials  and  alloys.  The  model  presented  here  is  a first  step  to  developing  a phase-field 
model  for  the  solidification  of  an  alloy.  In  future  work  we  hope  to  incorporate  this  model 
in  a formulation  based  on  increcising  entropy,  which  will  then  allow  both  variations  in  tem- 
perature and  composition,  and  provide  a phase-field  model  for  commercially  important  alloy 
sohdification  processes  used  in  such  things  as  the  fabrication  of  modern  electronic  devices 
and  directionally  solidified  turbine  alloys. 

In  section  2 we  discuss  the  phase-field  model  of  a pure  material  due  to  Kobayashi  [16], 
and  use  it  as  a basis  to  derive  our  isothermal  phase-field  model  for  growth  of  binary  alloys. 
We  employ  a kinetic  formulation  which  ensures  that  the  Helmholtz  free  energy  decreases 
monotonically  in  time.  In  the  following  section  we  consider  the  sharp  interface  limit.  This 
shows  the  correspondence  of  our  phase-field  model  with  classical  sharp  interface  models,  and 
allows  us  to  relate  the  materials  parameters  in  the  phase-field  model  to  those  of  the  classical 
Stefan  problem. 

We  show  that  the  time  evolution  of  the  system  occurs  on  three  scales.  There  is  an  initial 
very  short  period  in  which  interfaces  form.  It  is  a response  of  the  system  to  the  initial  data 
and  does  not  necessarily  represent  the  process  of  nucleation.  This  is  followed  by  a second 
period  in  which  the  interfaces  move  in  response  to  the  energy  density  difference  across  them. 
In  the  third  and  final  period  coarsening  occurs,  in  which  the  interfaces  move  very  slowly, 
primarily  in  response  to  changes  in  curvature.  In  section  4 we  discuss  the  motion  of  a planar 
interface  and  show  how  segregation  is  dependent  on  the  interfacial  growth  rate  in  this  model. 
In  the  last  section,  we  present  results  from  a numerical  integration  of  the  model,  which  bear 
out  the  results  of  section  3 and  4 and  show  the  power  of  this  formulation. 
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2 Problem  Formulation 


2.1  The  Governing  Equations 


We  consider  an  isothermal  solution  of  two  different  components  A and  B which  may  exist 
as  two  distinct  phases,  solid  and  liquid,  contained  in  a fixed  region  Q with  boundary  dO.. 
We  denote  the  concentration  of  B by  c(x,  t)  and  we  introduce  a phase  field  (^(x,  t)  which 
represents  the  phase  in  time  and  space  in  Specifically,  we  describe  the  solid-hquid  interface 
by  <^(x,  t)  = 1/2  and  denote  regions  where  <^(x,  t)  > 1/2  as  sohd  and  regions  where  <^(x,t)  < 
1/2  as  Hquid. 

A recent  phase-field  model  due  to  Kobayashi  [16]  models  the  phase  transition  of  a pure 
material  by  employing  the  following  gradient  weighted  free-energy  functional: 


dn. 


(6) 


where  e is  a constant,  T(x,  t),  is  the  temperature  and  the  free  energy  density  is 

fi^,T)=f%{p-l){p-^-/3iT))dp.  (7) 

Jo 

Here  0{T)  is  a monotonic  increasing  function  of  T,  such  that  P{Tm)  = 0,  where  Tm  is  the 
freezing  temperature  of  the  material,  and  \P[T)\  < 1/2.  The  free  energy  density  f((t),T) 
is  a double-well  potential  in  (j>^  see  Figure  1.  The  restriction  \^{T)\  <1/2  ensures  that  it 
has  local  minima  at  <^  = 0 and  (^  = 1,  and  a local  maximum  at  = 1/2  Because 

of  the  two  minima,  the  system  may  exist  stably  in  a state  which  is  a single  phase  liquid 
((^(x,  i)  = 0)  or  a single  phase  solid  (<^(x,t)  = 1).  There  is  an  energy  penalty  for  a change 
of  phase  within  the  region  fl,  which  corresponds  to  (j)  varying  between  zero  and  unity.  This 
is  because  such  a variation  increases  the  total  energy  of  the  system,  due  to  an  increased 
energy  density  associated  with  the  double- well  nature  of  the  energy  density,  and  also  due  to 
the  contribution  to  the  total  energy  due  to  the  gradient  energy,  which  is  no  longer  zero. 

If  — 1/2  < /3  < 0,  i.e  T < Tmi  then  the  global  minimum  of  the  energy  density  is  at  (^  = 1, 
and  so  the  single-phase  sohd  has  the  lowest  energy.  However,  if  0 < /3  < 1/2,  i.e  T > Tm 
then  the  situation  is  reversed,  and  the  single-phase  hquid  has  the  lowest  energy. 

We  now  employ  this  form  for  the  free  energy  density  to  develop  the  appropriate  free 
energy  density  for  an  isothermal  solution  in  which  the  two  different  phases  are  assumed  to  be 


6 


ideal  solutions.  We  assume  that  where  and  are  the  freezing  temperatures 

of  pure  A and  pure  B when  c = 0 and  c = 1 respectively.  We  assume  that  the  actual 
temperature  of  the  solution,  T,  which  is  given,  hes  between  these  two  freezing  temperatures, 
i.e.  > T > T^.  We  also  assume  that  the  Helmholtz  free  energy  densities  of  each 

component  A and  B aJone  are  of  the  form  given  by  Kobayashi,  and  are  denoted  by  /^((^;  T) 
and  respectively.  SpecificaJly  we  put: 

fA{cl>;T)  - Wa  [%{p  - l){p  - i - fiA{T))dp,  (8) 

Jo 

fB{<l>;T)  = Wb  f%{p  - l){p  - - MT))dp,  (9) 


where  Wa,Wb  are  constants,  ajid  the  temperature  T is  a parameter  in  this  isothermal 
situation.  We  note  that  because  < T < then  —1/2  < Pa{T)  < 0 < Pb{T)  < 1/2. 
We  take  the  energy  density  /(<^,  c;T)  of  the  solution  as: 


f{4>,c-,T)  = c/b(i^;T)  + (1  - c)fA{^;T)  + [c  logc  + (l  - c)  log(l  -c)],  (10) 

I'm 

where  R is  the  universal  gas  constant  constant  and  Vm  is  the  molar  volume,  which  is  assumed 
to  be  constant.  The  first  two  terms  correspond  to  the  contribution  to  the  energy  density 
due  to  the  individual  molar  Helmholtz  free  energies  densities  of  the  two  components  and  the 
last  term  is  due  to  the  decrease  in  energy  associated  with  the  mixing  of  the  two  components, 
under  our  assumption  that  both  liquid  and  sohd  phases  are  ideal  solutions.  In  Figure  2 we 
represent  the  free  energy  density  (10)  as  a surface.  This  model  could  easily  be  extended  to 
the  case  of  a regular  solution  by  adding  a term  of  the  form  c(l  — c)[n£,  + — Ql)]  to  the 

free-energy  density  /(c,  ^),  where  and  Qs  are  thermodynamic  constants. 

We  define  the  free  energy  functional  by: 


n^,c;T]=  f 

Jn 


dCl. 


(11) 


In  order  to  derive  a kinetic  model  we  make  the  assumption  that  the  system  evolves  in  time 
so  that  its  total  free  energy  decreases  monotonically.  The  following  prescription  is  consistent 
with  this  cLssumption  and  ensures  that  the  total  amount  of  solute  is  conserved: 


d<j> 

It 


(12) 


7 


| = A/.V.(<l-c)vf 


(13) 


where  Mi  and  M2  are  positive  constants.  The  choice  of  boundary  conditions 


d<l>  _ dc 
dn  dn 


(14) 


where  n is  the  outward  normal  to  the  boundary  50,  ensures  that  there  is  no  change  in  the 
total  composition  of  the  system  due  to  transport  across  its  boundary  and  that  the  Helmholtz 
free  energy  decrecises  monotonicaJly  with  time.  We  may  interpret  the  right  hand  side  of  (13) 
as  the  negative  of  the  divergence  of  a solute  flux,  j = — M2c(l  — c)V{8!F / 5c).  The  coefficient 
c(l  — c)  has  been  included  to  ensure  that  the  diffusion  equation  for  the  solute  that  emerges 
has  a diffusion  coefficient  that  is  not  dependent  upon  composition. 

Evaluating  the  variational  derivatives  of  the  free  energy  functional  gives: 


dt 

dc  _ 
dt 


"■V 


) 


which  may  be  also  written  by  using  (10)  as: 


dt 


= Mi 


dc 

dt 


M2V  • (c(l 


d<t> 


d4> 


c)V{fA  - fa))  + DV^C, 


(15) 

(16) 


(17) 

(18) 


where  D = M^RTIvm  is  identified  as  the  diffusivity  of  the  solute.  We  note  that  in  this 
formulation  the  solute  diffusivity  is  the  same  in  each  phase;  this  restriction  could  be  overcome 
by  allowing  the  quantity  M2  to  depend  on  <!>. 

We  now  nondimensionalise  these  governing  equations  by  introducing  the  following  di- 
mensionless quantities  (denoted  by  a tilde): 

7*2  nm 

x = rx,  t=—i,  T = T*f,  / = — /,  (19) 

IJ  Vm 

where  I*  is  a representative  length  of  the  system  12,  T*  = {T^  -}-  3j§)/2  is  the  average  of  the 
melting  temperatures  of  the  two  constituents,  and  time  has  been  scaled  on  the  diffusive  time 
scale.  Here  / = /(<^,c;  f)  = cfB{(j>\f)  + (1  - c)fA{^;  f)  -f  [c  logc-f-  (1  - c)  log(l  - c)],  where 
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fAi4>;  f)  = Wa  So  pip  - l)(p  - 5 - PA{f))dp,  fsi^’,  f)  = Wb  So  p{p  “ 1)(p  “ 5 - PB{f))dp, 
Wa  = Wav^IRT  and  Wb  = Wbv^IRT. 

The  governing  equations,  (17)  and  (18)  can  then  be  written  as 


where 


Ml  = 


DVrrv 


. _ £ ! Vm 


I*  V RT* 


(20) 

(21) 

(22) 


2.2  Relationship  Between  Materials  Parameters  of  the  Phase- 
Field  Model  and  Sharp  Interface  Theories 

In  order  to  determine  the  relative  magnitudes  of  the  dimensionless  parameters  and  to  per- 
mit application  to  real  materials  systems,  it  is  necessary  to  obtain  estimates  for  the  seven 
materials-related  parameters  of  the  alloy  phase-field  model.  They  are  e,  Wa,  Wb,/3a{T),/3b{T), 
Ml  and  M2.  The  parameter  M2  has  already  been  related  to  the  diffusion  coefficient  as 
D = M2RT jvm,  through  the  similarity  of  (18)  (when  (/>  is  identically  zero  or  unity)  to  the 
linear  diffusion  equation. 

We  proceed  to  develop  expressions  that  relate  the  parameters  e^Wa, Pa{T)  and  Mi  to 
parameters  common  to  sharp  interface  theories  for  pure  A.  These  parameters  are:  a a,  solid- 
liquid  surface  energy  of  a stationary  planar  interface;  La,  the  latent  heat  of  fusion  per  unit 
volume;  and  the  linear  kinetic  coefficient,  which  relates  the  normal  interface  velocity 
to  interface  temperature  by 

K = HAiT^  - T).  (23) 

This  data  (or  estimates  of  it)  is  available  for  most  pure  materials.  The  corresponding  iden- 
tification of  parameters  for  the  case  of  pure  B is  completely  analogous.  We  consider  each 
parameter  separately  in  turn. 
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2.2.1  Energy  Difference  Between  Phases 


The  difference  in  free  energy  between  solid  and  liquid  phases  is  known  as  a function  of 
temperature.  This  quantity,  /^(l;  T)  — /a(0;  T),  is  often  approximated  by  a linear  function 
of  temperature  near  the  melting  point  as: 


h{l;  T)  - Ua-,  T)  = 

Evaluation  of  this  free  energy  difference  from  the  integration  of  (8)  gives 

WaPa{T) 


(24) 


hence 


fAil-,T)-fA{0-,T)  = 


WaMT)  La{T-T^) 
6 


(25) 


(26) 


2.2.2  Free  Surface  Energy 

An  expression  for  the  excess  free  energy  of  the  interfacial  region  over  the  bulk  phases 
has  been  derived  for  an  interface  with  a nonconserved  order  parameter  (the  case  employed 
in  our  phase-field  model)  by  Allen  and  Cahn  [19].  With  change  of  notation  this  is 


= e 'J fA{<i>-,T)d4>. 


(27) 


Evaluation  of  this  integral  (at  T = Tm)  using  the  form  of  the  free  energy  density  given  by 
(8)  yields 


cja  = 


12 


(28) 


2.2.3  Interface  Kinetics 


We  consider  the  one- dimensional  solution  of  (15)  for  a pure  material,  c = 0 or  c = 1,  (see 
Hajowell  and  Oxtoby  [20])  which  represents  the  uniform  motion  of  a planar,  diffuse  interface. 
Transforming  to  a frame  moving  with  constant  velocity  Ki,  (15)  has  a steady  solution  which 
satisfies 

dx^  Ml  dx  d<l> 
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The  only  solution  with  the  property  that  ^^lascc— >— oo  and  0 0 as  cc  — ^ +oo  is 

T -1 


^(i)  = 


l + exp\xJ^ 


(30) 


and  this  solution  only  exists  when 


V„  = -Mj6/3A(T)y/^. 


(31) 


We  see  from  the  form  of  the  solution  (30),  or  from  Allen  and  Cahn  [19],  that  the  thickness 
of  the  interface  Sa  is  characterised  as 


Sa  = 

Equations  (26)  and  (31)  give 

Vn  = -6MieLA‘ 


Wa 


(32) 


W, 


qnA 


(33) 


from  which  by  comparison  with  (23)  we  find 

6MiLaSa 


= 


ta 

-L  HA 


(34) 


M 


Equations  (26),  (28)  and  (34),  derived  from  comparison  with  classical  sharp  interface 
theory,  provide  three  relations  for  the  four  unknowns  e,  W^,/?>i(r)  and  Mi.  In  addition, 
information  is  required  about  either  the  height  of  the  maximum  of  the  double- well  potential 
(I40i/64  at  T = T^)  or  the  interface  thickness  8a  to  uniquely  establish  values  for  all  pa- 
rameters. In  some  cases,  such  as  spinodal  decomposition,  the  height  of  the  maximum  of  the 
double-well  potential  can  be  estimated  from  a regular  solution  model  of  the  phase  and  the 
proximity  of  the  critical  point.  In  the  present  application  to  liquid- solid  reactions,  such  an 
approach  is  not  currently  possible.  We  therefore  choose  to  treat  the  interface  thickness  as 
an  input  parameter,  in  which  case  the  values  of  the  phase  field  are  related  to  the  classical 
data  as  I 

(35) 


Wa  = 12^2°^^ 

JWi  = 


Sa 


6La5a’ 


(36) 

(37) 
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(38) 


( LaSa  \ (T-T^\ 

\2v^(jyi/  V ) 

by  employing  equations  (26),  (28),  (32)  and  (34). 

We  emphasize  that  whatever  choice  we  mcike  for  the  interface  thickness  the  resulting 
values  of  e,  Wa,/3a(T')  and  Mi  provide  a model  that  has  the  required  free  energy  change, 
surface  energy  and  kinetics  to  match  the  classical  theory. 

To  treat  the  alloy  case  at  hand  in  the  present  paper,  the  last  four  equations  would  be 
evaluated  for  each  pure  component.  However,  because  we  have  assumed  that  e and  Mi 
are  constant  and  do  not  depend  on  composition,  we  are  restricted  in  the  choices  for  some 
parameters.  To  ensure  that  single  values  of  e and  Mi  are  obtained,  the  following  restrictions 
must  be  satisfied: 


— — — 
o-A  ’ 

fJ'A  _ La 
fiB  ^B  Lb 


(39) 

(40) 


These  restrictions  could  be  removed  by  choosing  both  e and  Mi  to  depend  linearly  on  com- 
position. 

Numerical  values  of  these  constants  can  be  estimated  for  the  alloy  system  Ni-Cu.  The 
liquid-solid  part  of  this  phase  diagram  closely  resembles  that  for  the  ideal  solution  model 
assumed  in  this  work.  Values  for  the  parameters  are  taken  primarily  from  Coriell  and 
Turnbull  [22]  and  are  given  in  Table  1.  The  dimensionless  parameters  e and  Mi  can  be 
calculated  from  the  data  in  Table  1.  For  a length  scale  /*,  we  chose  a dimension  typical  of 
cellular  and  dendritic  structures  of  100/zm.  Using  the  data  in  Table  1 and  T*  = (2m  + Tj^)/2, 
we  have  e = 9.6  x 10~®,  and  Mi  = 8.5  x 10^^.  These  values,  along  with  equations  (20) 
and  (21),  provide  a complete  formulation  for  the  calculation  of  isothermal  phase  transition 
(solidification)  in  the  Ni-Cu  alloy. 


3 Asymptotic  Analysis:  The  Limit  e — 0. 

In  this  section  we  analyse  the  system  in  the  limit  e — ► 0.  We  restrict  our  attention  to  the 
situation  when  the  domain  Q is  two-dimensional.  We  assume  that  (j)  and  c satisfy  Neumann 
conditions  on  the  boundary,  so  that  the  total  amount  of  solute  is  conserved  and  the  free 
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energy  of  the  system  decreases  monotonically  in  time.  In  particular,  we  consider  the  spatio- 
temporal  evolution  of  the  concentration  and  phase  field  from  given  initial  data.  We  require 
that  the  initial  data,  Ctntt(x)  and  for  c(x,  t)  and  ^(x,  t)  respectively  is  bounded  in  the 

interval  [0,1],  i.e.  0 < Cinit(^)  ^ Ij  0 < <^mtt(x)  < 1,  and  satisfies  the  Neumann  boundary 
conditions.  We  emphasize  that  because  of  the  generality  of  this  initial  data  it  does  not 
necessarily  represent  a physically  realistic  initial  state. 

We  assume  that  the  contribution  of  the  gradient  energy  is  much  smaller  than  that  of 
the  free  energy  density  and  so  « RT*  jvyn^  and  hence  e <<  1.  Further,  we  assume 

Ml  >>  1 and  in  particular  we  put: 

M 

Ml  = — , (41) 

e 


where 

D \ Vm 

In  fight  of  the  values  of  Mi  and  e calculated  above  for  the  Ni-Cu  alloy  such  assumptions  are 
clecirly  justified. 

The  governing  equation  for  the  phase  field  becomes: 


(43) 


and  that  for  the  composition  is  unchanged  from  (21).  Below  we  shall  consider  the  limit 
e — ^ 0,  which  corresponds  to  considering  the  situation  when  the  mobility  of  the  phase  field  is 
very  large  and  the  thickness  of  the  interfaces  is  very  thin.  We  note  that  this  represents  one 
distinguished  limit  of  the  governing  equations  (21)  and  (43),  which  we  show  below  provides 
a realistic  model  of  phase  transitions  in  isothermal  alloys.  Other  distinguished  limits  are  of 
course  possible,  as  for  instance  considered  by  Caginalp  [12],  in  his  analysis  of  the  phase-field 
model  of  a pure  material. 

In  the  limit  e ^ 0 the  diffusivity  associated  with  the  phase  field,  which  is  Me,  becomes 
small,  cind  results  in  the  formation  of  interior  layers  within  the  region  Q,  in  which  the  phase 
field  and  the  concentration  vary  rapidly,  the  former  from  zero  to  one.  Such  layers  represent 
solid-liquid  interfaces.  In  this  section  we  seek  to  derive  the  jump  conditions  which  relate  the 
concentration  field  and  its  derivatives  in  the  vicinity  of  such  layers  to  their  normal  velocity 
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and  curvature.  We  first  consider  the  outer  solution  in  regions  away  from  the  interfacial 
layers. 

Henceforth  we  only  deal  with  dimensionless  quantities  in  the  governing  equations,  and 
for  simphcity  of  notation,  we  omit  the  tilde’s. 

3.1  Outer  Solution 

To  proceed  we  seek  an  outer  solution  as  a regular  perturbation  in  e.  Thus  we  put: 


^(x,<)  = + e^(^)(x,f)  + C>(e^),  (44) 

c(x,t)  = + 0{e^).  (45) 

The  leading  order  problem  is 

|f(^W.c(«))  = 0.  (46) 

— = Vc(°)(l  - c(°))  • V(/^  - /b)  + (47) 


from  (46)  we  find  that  = 0 or  = 1,  depending  on  whether  the  region  is  all  Hquid 
or  all  solid  respectively.  (The  other  extremum  of  / represents  an  unstable  maximum  of  the 
energy,  and  need  not  be  considered.)  It  follows  then  from  (47)  that  the  leading  order  solute 
field  simply  satisfies  the  diffusion  equation: 

^ (48) 

The  governing  equations  for  the  first  order  corrections  may  also  found  in  a similar  way; 
however,  it  is  easily  shown  that  there  are  no  contributions  to  <l>  beyond  the  leading  order 
terms,  i.e.  is  exactly  zero  or  one  as  appropriate. 

Because  of  the  generahty  of  the  initial  data  that  we  have  allowed,  the  initial  conditions 
will  not  necessarily  correspond  to  (j)  zero  or  unity.  In  such  a situation  it  is  simple  to  show 
that  there  is  an  initial  fast  transient,  on  a time  scale  0(e),  in  which  <!>  attains  the  value 
zero  or  unity  for  cdmost  each  x E depending  on  the  value  of  the  initial  data  for  (f)  at  x. 
After  this  short  period  there  will,  in  general,  be  large  gradients  in  <j),  which  represent  the 
interfacial  layers.  During  this  short  initial  period  interfaces  are  created.  It  is  their  subsequent 
development  which  we  now  consider. 
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3.2  Inner  Solution:  The  Interfacial  Layer 


We  assume  that  the  interfacial  layer  ((^(x,t)  = 1/2)  is  given  by  the  curve  r(t).  We  introduce  a 
local  orthogonal  coordinate  system  (^(x,  t),  s(x,t))  where,  in  the  vicinity  of  r(t),  r measures 
the  normal  distance  to  the  curve.  Thus  T{t)  is  therefore  described  by  r(x,t)  = 0,  and 
|Vr|  = 1 in  its  neighborhood.  Further  we  note  that  the  curvature  )C  of  T[t),  is  given 
hy  )C  = V^r.  Along  r(t)  the  coordinate  ^(x,  t)  measures  distance  from  some  fixed  point. 
Employing  these  new  independent  variables  the  governing  equations  become: 


(t>t  + rt(t)r  + St<f>3  = M e[(l)rr  + - € 


Ct  + TtCr  + 5tC,  = 


+c(l  - c))  j 

(tc 


,2,df 


I) 

/ rr  \ y 33 


+ W 


d<f> 


dc 


(49) 


(50) 


In  the  interfacial  layer  the  rapid  changes  in  the  phase-field  variable  ensure  that  the  effect 
of  the  diffusion  is  important.  Thus  we  scale  the  coordinate  r with  e as: 


r = ep. 


We  also  introduce  an  additional  slow  time  scale,  r: 


T = ct. 


(51) 


(52) 


We  shall  employ  the  method  of  multiple  scales  and  hence  we  associate  the  partial  derivative 
with  respect  to  time  in  the  governing  equations  with  d/dt  edj dr. 

We  expand  the  dependent  variables  as: 

c=  ^°\p,3,t,T)  + S^\p,s,t,T)  + 0{e^),  4i  = $^°\p,3,t,T)  + e^^^\p,  s,t,T)  + (53) 

where  we  denote  variables  in  this  region  by  an  over  bar. 

The  matching  conditions  upon  the  inner  solution  in  the  transition  region  are  found  by 
the  expanding  the  outer  solution  in  terms  of  the  inner  variables, 

ac(°) 


C — ^ C^°^|(r=0+)  + C + P ■ 


+ 0{e^),  as  p ^ CO, 


(54) 


•=o+. 
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c ->  + e ( c(^^|r=o-  + 


dr 


+ O(e^),  as  /?  — > — oo, 


(r=o-), 

^ — 1 duS  p — OO, 


(^  = 0 as  ^ +00. 


(55) 

(56) 

(57) 


In  choosing  the  matching  conditions  in  this  way  we  have  assumed  that  in  this  local  coordinate 
system  the  outer  liquid  region  is  aJiead,  and  the  outer  solid  region  is  behind,  the  interfacial 
region.  To  proceed  we  insert  the  forms  (51),  (52)  and  (53)  into  the  governing  equations 
(49)  and  (50)  to  yield  a sequence  of  partial  differential  equations  at  each  order  of  e,  which 
we  now  consider  in  turn. 


3.2.1  Leading  Order 
The  leading  order  problem  is: 

= M ^(°)))  . 

=0. 

The  latter  may  be  integrated  to  give  that: 

c^°^(l  — c^°l)  j — - = function  of  s,  t and  r, 

dp 

However  the  matching  conditions  (54),  (55),  (56),  (57),  require  that  dfc/dp 
Hence  we  deduce  that  the  function  of  integration  is  zero  and  so: 

fc{^°\$^°^)-A{s,t,T)  = 0,  (61) 

where  A is  a function  of  integration,  which  at  this  stage  is  undetermined.  It  is  now  convenient 
to  define  a new  modified  energy  density  ^^°^)  by 

^(“))  = /(c^o),  ,^(“))  - A{s,  t,  t)c(°).  (62) 

in  which  case  the  equations  (58),  (59)  may  be  expressed  as: 

= M (0(°)  - F*(c<“),  ^(“)))  (63) 


(58) 

(59) 

(60) 

0 as  \p\  oo. 


16 


F,(c(°),^W)  = 0. 


(64) 


By  employing  the  matching  conditions  we  find 

Fcic^°\0-,t),l)  = Fc(c‘°>(0+.4).0)  = 0.  (65) 

Equation  (64)  relates  to  as: 

= ffo(0'“’;  A),  (66) 

where, 

So(^^°';  A)  = (l  + exp  {/s(i^(°>)  - /x(^'°*)  - A(5,  i,  r)})  \ (67) 

Thus,  the  modified  energy  density  F(c(°l,  A),  is  a function  of  0^°) 

only  which  we  denote  by  (=  ^)j  in  which  case  (63)  and  (64)  may 

be  combined  as  the  single  equation: 

« - ^ (“) 

This  problem  determines  (and  hence  c^°l  from  (66))  and  i4(s,  t,  r).  The  function  A(s,  t,  r) 
is  found  from  the  asymptotic  behaviour  required  by  the  matching  conditions,  which  is  only 
possible  if  the  boundary  conditions  represent  singular  points  of  (68),  in  which  case  the 
required  solution  is  the  heteroclinic  orbit  between  them.  Any  other  solution  would  attain 
the  boundary  conditions  at  a finite  value  of  p. 

We  now  derive  an  expression  for  the  normal  velocity  Vn  = —rt.  We  multiply  (68)  by 
and  integrate  to  obtain: 

^ /_  = 0-  (69) 

However  from  matching: 

^<,°>U=o  = 4°>|p=-cc  = 0,  (70) 

and 

» 0 as  p — > oo  (71) 

— > 1 as  p — oo  . (72) 
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and  so: 


M[E{0)  - £;(!)] 


(73) 


where 

To  = r (74) 

J — OO 

Thus  we  see  that  the  interface  motion  is  driven  by  the  energy  disequilibrium  between  the  two 
phases;  the  interface  locaJly  advances  into  the  phase  with  the  higher  modified  energy  density, 
The  interface  acts  to  reduce  the  energy  F of  the  system  by  moving  to  convert 
regions  of  high  F into  low  F.  Time  only  enters  the  problem  for  through  the  matching 
conditions  and  hence  acts  as  a parameter.  Thus  and  Vn  are  time  dependent  and  are 

controlled  by  the  time-dependence  of  the  solute  field  outside  the  interfacial  layer.  In  section 
4 we  show  that  equations  (65)  and  (73)  may  be  interpreted  geometrically  as  a parallel 
tangent  construction  to  the  Helmholtz  free-energy  surface  for  a moving  planar  interface  or 
a common  tangent  construction  for  a planar  stationary  interface. 


3.2.2  First  Order 


We  now  go  on  to  consider  the  next  order  approximation  to  the  governing  equations: 

(75) 

= (c(“>(l  - ^°')(/cc(c<“>,  . (76) 

We  may  integrate  (76)  once  from  ^ = — ooto  /o  = +oo,  which  gives: 


p=oo 


p=oo 


= 0. 


(77) 


\ dr 


\ = ~Vn  (c(°^|(r=o+)  - c(°^|(r=0-))  • 
0-/ 


(78) 


On  employing  the  matching  conditions  this  gives  that: 

dr 

0+  ^ 

This  simply  represents  a conservation  of  solute  normal  to  the  moving  interface  and  is  the 
usual  condition  employed  in  the  classical  model  model  of  phase  transitions  involving  sharp 
interfaces. 
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As  an  isolated  system  evolves,  it  approaches  its  final  minimum  energy,  in  which  case 
the  energy  differences  across  the  interface,  which  drive  the  interfacial  motion,  become  small. 
Thus  > 0 as  t — > CO  and  so  the  system  evolves  on  a slower  time  scale.  To  determine 
the  interfaciaJ  d5mamics  on  this  long  time  scale,  which  is  given  by  the  scaled  time  r,  we 
put  T4  = 0,  in  which  case  the  interfacial  velocity  is  given  by  In  this  regime  the 

concentration  gradients  near  the  interfaces  are  small,  0{e).  The  slow  velocity  Vn  is  given 
from  a solvability  condition  that  arises  from  the  first  order  problem.  In  the  Appendix  we 
discuss  the  adjoint  problem  and  use  it  to  determine  the  solvability  condition,  which  is: 


c(o)(0+)-c(°)(0-) 


+ ro(Af-v„  + x:)  = 0, 


(79) 


c(o)(0-)(l-c(o)(0+)) 

where  Fq  is  given  by  equation  (74).  The  above  solvability  condition  determines  the  slow 
component  of  the  normal  interface  velocity,  which  is  thus  dependent  upon  the  local 
curvature,  /C,  and  the  perturbed  concentration  in  the  Hquid,  In  the  Appendix 

we  also  derive  another  solvabihty  condition,  equation  (116),  which  relates  the  perturbed 
interfacial  concentrations  c(^)(0“)  and  c(^)(0''‘)  at  the  interface. 


3.2.3  Temporal  Evolution 

From  the  above  analysis  we  can  identify  three  different  stages  in  the  evolution  of  the  solution 
of  the  phase-field  equations  based  on  initial  data  in  which  0 < (j)init  £ 1: 

• Stage  I.  Initially  there  is  a very  fast  initial  transient  0(e)  in  which  the  phase  field 
adjusts  locally  to  the  values  zero  or  unity.  We  interpret  this  as  a period  in  which 
interfacial  layers  are  born  and  regions  of  solid  and  liquid  are  differentiated. 

• Stage  II.  Subsequently,  on  an  0(1)  time  scale,  away  from  these  interfacial  layers  in 

the  outer  regions,  (f)  is  either  zero  or  unity  and  the  solute  concentration  is  governed  by 
the  diffusion  equation.  Between  the  outer  regions  the  interfacial  layers  exist.  The  local 
normal  velocity  of  these  interfaces  is  determined  by  the  local  jump  in  the  modified 
energy  density  across  the  layer,  E,  (and  thus  by  the  deviation  of  the  composition  from 
the  equilibrium  values)  described  by  (73).  More  specifically,  the  interfaces  move  so 
as  to  decrease  the  energy  of  the  system.  The  interfaciaf  motion  will  slow  down  as  the 
system  approaches  its  minimum  energy,  i.e.  — > 0 as  t > oo. 
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• stage  III.  Finally,  on  a longer  time  scale  0(e~^),  the  normal  interfacial  velocity  is 
given  from  (79).  The  slow  motion  of  the  interface  is  controlled  by  a balance  between 
curvature  and  thermodynamic  disequilibrium,  the  latter  is  represented  by  the  term 
involving  the  perturbed  concentration.  This  represents  the  process  of  Ostwald  ripening. 

Thus,  initially  the  system  evolves  very  fast  to  minimise  its  energy,  by  evolving  the  phase 
field  at  each  point  in  space  to  zero  or  unity;  the  minima  of  the  energy  density  for  a fixed 
concentration.  Subsequently  diffusion  acts  and  the  interfaces  move  to  reduce  the  energy 
differences  across  them.  Finally  they  slowly  evolve  in  response  to  local  energy  disequilibrium 
and  kinetic  effects,  as  well  as  to  reduce  the  energy  associated  with  their  curvature. 

3.2.4  A Mechanical  Analogy 

The  leading  order  behaviour,  discussed  above  is  described  by  the  system: 

+ (80) 

0 = F^{x,y),  (81) 

when  we  identify  x with  y with  p with  time,  and  F{x,y)  = f{x,y)  — Ax  as 
before,  where  A is  to  be  determined  and  / is  the  free  energy  density.  This  system  represents 
the  motion  of  a paxticle  in  a potential  z = —F{x,y),  constrained  to  lie  on  the  curve  x = 
go{y]  A),  with  friction  acting  only  in  the  y direction.  We  note  that  the  surface  z = F(x,y)  is 
obtained  from  that  of  the  free  energy  density  surface  z = /(x,t/)  by  inverting  it  and  tilting 
it  in  the  x direction;  the  amount  of  tilt  depending  on  the  value  of  A.  The  frictional 

constant  is  equal  to  M~^Vnj  and  so  the  stationary  interface  corresponds  to  the  frictionless 
case.  The  matching  conditions  require  that  the  particle  follows  a path  such  that  it  approaches 
a stationary  state  as  time  goes  to  plus  and  minus  infinity.  The  only  possible  such  motion  is 
one  that  begins  and  ends  at  local  maxima  of  the  surface  z = — F'(x,t/).  The  form  of  the  free 
energy  density  ensures  that  such  local  maxima  only  occur  for  one  given  value  of  x on 

y = 0 = 0)  and  one  given  value  of  x on  y = I = 1),  both  of  which  depend 

on  A and  He  on  the  curve  x = go{y]  A)  . Thus  the  particle  must  follow  a trajectory  between 
these  local  maxima  which  traverses  from  y = 0 to  y = 1.  When  the  interface  is  stationary. 
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the  particle  moves  in  the  absence  of  friction,  in  which  case  the  constant  A (or  equivalently 
the  amount  of  tilt  of  the  surface)  is  chosen  such  that  the  heights  of  the  two  maxima  are  the 
same.  For  a non- station  ary  interface  the  constant  A must  be  chosen  such  that  the  potential 
energy  difference  associated  with  the  difference  in  height  of  the  two  maxima  is  just  sufficient 
to  account  for  the  energy  dissipated  by  friction.  A mechanical  analogy  was  noted  by  Cahn 
[4]  for  the  Cairn- Hilliard  equation,  as  well  as  by  Harrowell  and  Oxtoby  [20]. 

3.3  Discussion  Of  Asymptotics 

As  shown  in  section  3.1  the  governing  equation  for  the  solute  reduces  to  the  diffusion  equation 
in  the  solid  and  hquid  regions  with  the  diffusivity  D given  by  B = MiRT* /vm.  Thus  the 
governing  equation  for  composition  given  here  is  the  same  as  those  of  a sharp  interface,  or 
classical  model,  in  the  solid  and  liquid  regions. 

We  see  from  (73)  that  the  temperature,  T,  and  normal  velocity,  Vn,  determine  the 
leading  order  interfacial  concentrations  in  the  solid  and  liquid.  Thus  cs  = HsiT^Vn)  and 
Cl  = HLiT^Vn).  When  the  interface  is  stationary,  = 0,  then  the  functions  cs  — Hs{T^^) 
and  Cl  = HL{T^{f)  simply  represent  the  solidus  and  liquidus  curves  of  the  phase  diagram  of 
the  alloy  (see  section  4.1). 

The  asymptotic  analysis  generalises  the  results  given  in  section  2.1  for  the  interface 
kinetics  and  surface  energy  of  a pure  material  to  the  case  of  an  alloy.  In  particular  (73)  is 
the  non-dimensional  analogue  of  (31)  for  the  kinetics  of  a pure  material.  The  solvability 
condition  (79)  represents  an  energy  balance  at  0{e)]  the  first  two  terms  are  the  change  in 
energy  due  to  the  perturbed  concentration  and  interfacial  motion,  the  last  term  represents 
the  surface  energy.  Thus  the  surface  energy,  a,  is  e JToo{i^p^Y^P  which  may  be  rewritten 
with  the  help  of  (69)  as 

a = eV2  / (82) 

which  is  the  generalisation  of  (27)  and  a similar  result  of  Cahn  and  Hilhard  [4].  However,  for 
the  alloy  these  expressions  for  the  surface  energy  (82)  and  interface  kinetics  (73)  have  no 
simple  analytic  form  as  they  do  for  the  case  of  the  pure  material.  For  simplicity  we  therefore 
choose  to  determine  the  parameters  of  the  phase-field  model  in  the  way  discussed  in  section 
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2.2  by  making  contact  with  the  properties  of  the  pure  components. 


4 Constant  Velocity  Planar  Solutions  as  e — > 0. 

In  this  section  we  look  for  solutions  in  the  limit  e — > 0 corresponding  to  a planar  interface 
advancing  at  at  constant  normal  velocity  V.  We  only  consider  the  outer  solution.  We  take 
a coordinate  coincident  with  the  moving  interface,  which  is  then  given  by  the  plane  z = 0. 
The  leading  order  problem  for  c is  then: 


Czz  + Vc^  = 0 
M[E{0)  - E{1)] 


To 


c|z=o~  c|^— 0+  ^o(0>  •^) 


dr 


z=0+ 


dc 

dr 


— —y{c\z=iO+  — c|j:=o-) 


(83) 

(84) 

(85) 

(86) 


z=0~  • 


with  0(2)  Coo  as  l^l  — ^ 00.  The  leading  order  solution  is  easily  found  to  be: 


c{z)  = C5  + (cl  — C5)exp(— Fz),  z>  0,  (87) 

c{z)  = cSiZ<0,  (88) 

where  we  associate  cs  and  cl  with  c(°^|z=o+  respectively,  which  along  with  the 

constant  A are  determined  by  the  three  conditions  given  in  (84)  and  (85).  We  note  that 
Coo  = Cs,  which  therefore  cannot  be  prescribed  independently  of  the  temperature,  T. 

The  segregation  coefficient  is 

jfe  = — , (89) 

Cl 

which  depends  on  the  interface  speed  V.  For  this  reason  we  first  consider  the  situation  when 
the  interface  is  stationary. 


4.1  Stationary  Interface:  V = 0 

In  this  case  the  concentration  is  constant  in  each  phase: 

c{z)  = cs,z  < 0,  c{z)  = cl,z>  0,  (90) 
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where  cs^cl,  and  A satisfy  the  three  conditions  in  (84)  and  (85)  with  V=0.  The  condition 
(84)  gives  that: 

i^(c5, 1)  = ^(cl,0).  (91) 

We  note  that  from  the  definition  of  the  function  go  given  by  (66)  and  the  conditions  in  (85) 
that: 

Fc(c5,1)  = Fc(cl,0)  = 0.  (92) 

Using  the  fact  that  F{cj(j>)  = — Ac  it  is  easily  shown  that  these  statements  are 

equivalent  to  requiring  that  cs  and  cl  are  chosen  such  there  is  a common  tangent  between 
the  curves  /(c,  0)  and  /(c,  1).  This  construction,  see  Figure  3(a),  is  precisely  that  given 
by  classical  thermodynamics  to  determine  the  equihbrium  concentrations  in  this  situation. 
It  is  usually  derived  by  requiring  that  the  chemical  potential  of  the  two  phases  on  either 
side  of  the  interface  in  a first  order  phase  transition  is  continuous.  The  locus  of  pairs  of 
compositions  cs  and  cl  at  different  temperatures  provides  the  hquidus  and  soHdus  curves 
of  the  phase  diagram.  The  common  tangent  construction  also  corresponds  to  constructing 
a common  tangent  plane  to  the  free  energy  surface,  shown  in  Figure  2,  which  is  tangent  to 
both  sections  of  the  energy  surfaces,  <^  = 0 and  <f)  = 1. 

We  now  employ  this  stationary  interface  solution  to  examine  the  possible  stationary 
configurations  of  the  system  as  e — > 0.  In  particular  we  consider  a one  dimensional  box, 
of  length  L,  whose  ends  are  impermeable  to  both  components  of  the  solution.  The  average 
composition  C is  ^ Jq  cdx.  There  are  single  phase  solutions;  single  phase  sohd  [(j)  = 1)  and 
single  phase  liquid  {(j)  = 0).  These  exist  for  aU  values  of  the  average  composition,  C G [0, 1], 
and  are  linearly  stable.  However,  the  analysis  above  indicates  that  other  solutions  are  possible 
representing  a single  interface  between  two  regions  of  different  phase  when  min  (c5,C£,)  < 
C < max(c5,  C£,).  In  Figure  4 we  schematically  illustrate  the  possible  solutions.  The 
variation  of  the  fraction  sohd  is  given  classically  (disregarding  the  interfacial  surface  energy) 
by  the  “lever”  rule  and  is  hnear  in  the  composition  between  cs  and  cl. 
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4.2  Non-Stationary  Interface  V ^ 0 

In  this  situation  the  only  condition  that  is  different  is  (91)  which  becomes: 

F(cl,  0)  - F{cs,  1)  = M-^VTo,  (93) 

the  left  hand  side  of  which  is  a constant,  AF{VyT),  for  a given  interface  velocity,  V,  and 
temperature,  T.  The  conditions  (92),  (93)  cire  equivalent  to  selecting  cs  and  cl  such  that 
their  tangents  on  the  curves  /(c,  0)  and  /(c,  1)  are  parallel,  but  the  tangents  are  a vertical  dis- 
tance AF( V,  T){=  F{cl^  0)  — F{cs,  1))  apart,  see  Figure  3(b).  It  follows  that  the  equilibrium 
concentrations  of  a moving  planar  interface  are  dependent  on  its  normal  velocity. 

From  (84)  and  (85)  we  find  that  the  segregation  coefficient  is  given  by: 

k = exp  {-AF{V,  T)  - fsil)  + /s(0)) , (94) 


which  for  small  V is  given  approximately  by: 

fc  = fco  (l  - ro|(v=o)  + <^(^'))  . (95) 


where  ko  is  the  value  of  the  segregation  coefficient  for  the  stationary  planar  interface.  Sim- 
ilarly for  a small  energy  difference  across  the  interface  the  segregation  coefficient  may  be 
approximated  by: 


k = ko 


1 - AF(0,  T)  -f  O((AF(0,  r))2) 


(96) 


From  equation  (68)  we  numerically  computed  c^,  cs  and  fc  as  a function  of  the  inter- 
face velocity.  The  results  are  displayed  in  Figure  5.  We  see  that  the  soHd-liquid  interfacial 
concentrations  and  the  segregation  coefficient  are  monotonically  decreasing  functions  of  the 
interface  velocity,  (the  latter  dependence  is  in  agreement  with  the  small  growth  rate  expres- 
sion for  fc,  (95)). 

Experimentally,  the  behaviour  of  the  segregation  coefficient  fc  with  increasing  interface 
velocity  has  been  shown  to  change  from  the  equilibrium  value  (normally  less  than  unity)  at 
low  velocities  to  values  necir  unity  for  high  velocities  ( « 1 m/s),  see  White  et  al.  [21].  This 
tendency  towards  unity  is  not  recovered  in  the  current  model.  We  anticipate  that  inclusion 
of  a term  in  the  energy  functional  involving  the  solute  gradient  energy  in  addition  to  the 
phase  field  gradient  energy  will  produce  a wide  variety  of  possible  behaviour,  depending  on 
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the  relative  importance  of  the  two  gradient  energy  terms.  In  fact  a model  by  Baker  [23] 
shows  that  k need  not  change  monotonically  towards  unity  with  increasing  velocity. 

5 Numerical  Calculations 

In  this  section  we  give  the  results  of  numerical  calculations  upon  our  model.  The  aim  of 
these  calculations  is  two-fold,  first,  to  verify  as  far  as  possible  the  asymptotic  analysis  given 
in  the  previous  section,  and  second  to  demonstrate  that  our  model  qualitatively  models 
the  behaviour  of  phase  transitions  in  an  isothermal  alloy.  In  order  to  apply  our  model 
to  situations  in  which  the  interface  is  sharp,  we  require  the  value  of  e to  be  small.  This 
provides  the  major  difficulty  in  conducting  calculations,  which  is  that  the  interfacial  layers  are 
necessarily  thin,  and  thus  a sufficiently  fine  computational  mesh  is  required  in  order  to  resolve 
the  solution  structure  within  them.  For  this  reason  it  is  much  easier  to  conduct  calculations 
in  one  spatial  dimension  than  two.  Thus,  in  this  section  we  first  employ  calculations  in 
one  dimension  to  verify  quantitatively  much  of  the  boundary  layer  analysis,  and  we  further 
perform  two-dimensional  calculations  with  a somewhat  larger  value  of  e to  provide  evidence 
that  our  model  is  qucditatively  in  agreement  with  this  analysis  and  that  it  shows  behaviour 
similar  to  that  found  in  real  systems. 

5.1  One-Dimensional  Calculations 

We  consider  the  evolution  of  the  system  (20)  and  (21)  in  one  spatial  dimension,  x,  on 
the  domain  0 < x < L.  We  imposed  the  boundary  conditions  = Cx  = 0 at  x = 0,  and 
<t>x  = 0,Cx  + aoc  = ai  a.t  X = L.  The  initial  data  is  given  by  (^(x,0)  = (^0(2:),  c(x,  0)  = co(x). 
We  used  the  software  package  PDECOL  [24],  which  employs  a finite  element  collocation 
procedure  for  the  spatial  discretisation  and  the  method  of  fines  to  advance  the  solution  in 
time.  In  all  cases  we  used  500  break  points.  The  first  set  of  calculations  were  aimed  to  test 
the  cLsymptotic  analysis  by  generating  the  planar  traveling  wave  solutions  predicted  by  the 
theory  given  in  section  4.  The  most  trivial  solution  corresponds  to  the  stationary  interface,  in 
which  case  the  asymptotic  analysis  predicts  the  (constant)  values  of  the  concentration  either 
side  of  the  interface.  With  initial  data  for  c and  (p  set  to  the  appropriate  step  functions. 
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it  is  found  that  the  predicted  concentrations  were  indeed  recovered  within  the  accuracy 
of  the  scheme,  with  e taken  sufficiently  small.  We  also  conducted  calculations  for  a non- 
stationary interface.  The  results  of  the  asymptotic  analysis  displayed  in  Figure  5 predict 
the  interface  velocity,  V to  be  0.1  when  cs  = 2.0339  x 10“^  and  C£,  = 0.3679.  In  order  to 
avoid  end  effects  we  chose  ao  = V and  ai  = Cg,  which  is  exactly  satisfied  by  the  traveling 
wave  solution  concentration  profile,  (87),  when  <^  = 0.  Physically  this  corresponds  to 
extracting  solute  from  the  liquid  phase  at  the  boundary  of  the  domain  in  order  to  maintain 
the  supersaturation  which  is  driving  the  interfacial  motion.  We  put  Pa  = —0.1, /3b  = 0.1, 
Wa  = Wb  = 1*0,  e = 10“^,  and  L = 0.4.  The  initial  data  was  <t)Q  = 0.1  + 0.9if(x  — 0.04) 
and  Co  = Cs  + H{x  — 0.04)(cl  — cs)exp[— F(x  — 0.04)],  where  H{x)  is  the  Heaviside  function. 
In  Figure  6 we  display  the  computed  composition  profile  at  times  of  0.0,  1.0,  2.0  and  3.0.  It 
clearly  reproduces  the  predicted  travehng  wave  solution  with  the  predicted  constant  velocity 
interface  velocity  of  0.1.  We  note  that  initially  the  velocity  of  the  interface  is  slightly  less 
than  that  predicted  by  the  asymptotic  analysis  as  e — > 0.  We  attribute  this  to  the  adjustment 
required  by  the  profile  in  the  interfacial  boundary  layer,  which  was  not  properly  given  by 
the  initial  data.  This  is  manifested  by  the  smoothing  out  of  the  initial  discontinuity  of  the 
initial  data  in  this  region.  The  above  calculations  quantitatively  verify  the  result  given  by 
(73)  and  confirm  that  the  modified  energy  differences  drive  the  interfacial  motion. 

We  also  conducted  a numerical  calculation  which  models  the  motion  of  a planar  interface 
in  a closed  system  in  response  to  the  presence  of  an  initially  supersaturated  liquid.  From 
Figure  5 the  interface  is  stationary  when  cs  = 0.1589  and  cl  = 0.8411.  We  chose  the  same 
parameters  as  the  previous  calculation  but  the  initial  data  for  the  composition  was  different: 
Co  = C5  + (0.6  — cs)H[x  — 0.04).  Also  the  boundary  condition  for  the  composition  at  x = 0.4 
was  Cx  = 0.  This  initial  state  corresponds  to  a supersaturated  liquid  phase  with  composition 
0.6(<  cb).  In  Figure  7 we  display  the  composition  profile  at  times  0,  2.0,  4.0,  6.0,  8.0  and 
10.0.  There  is  in  fact  a rapid  initial  transient,  not  apparent  in  this  figure,  in  which  the 
concentration  in  the  solid  decreases  to  approximately  0.1  and  the  concentration  in  the  liquid 
adjacent  to  the  interface  increases.  Subsequently,  as  shown  in  Figure  7 the  concentration 
boundary  layer  so  formed  in  the  hquid  is  accompanied  by  forward  interfacial  motion.  Finally, 
the  system  relaxes  back  to  the  uniform  concentrations  predicted  by  the  asymptotic  theory 
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for  a stationary  interface.  This  behaviour  is  expected  on  physical  grounds;  the  interface 
moves  to  reheve  the  supersaturation  subject  to  the  constraint  that  the  composition  of  the 
system  is  conserved.  However,  in  a real  system  it  would  be  expected  that  the  planar  interface 
would  be  unstable  and  evolve  to  a cellular  structure  as  it  moved  forward.  We  next  consider 
two-dimensional  calculations. 

5.2  Two-Dimensional  Calculations 

As  discussed  above,  the  thin  interfacial  layers  presently  cause  a difficulty  in  conducting 
accurate  numerical  calculations.  In  two  dimensions  this  problem  is  important  because  the 
computer  storage  and  processor  requirements  associated  with  fine  meshes  become  restrictive. 
Below  we  present  the  results  of  prehminary  calculations  which  employ  a uniform  80  x 80 
mesh,  the  finest  mesh  that  was  feasible  using  the  computer  facilities  available  to  us.  We  used 
the  software  package  PDETWO.  This  employs  second  order  accurate  central  finite  differences 
to  conduct  the  spatial  discretisation,  allied  to  an  accurate  ordinary  differential  integrator  to 
evolve  the  solution  in  time.  We  have  performed  a number  of  calculations.  Below  we  show 
typical  results  for  one  particular  case. 

We  performed  a calculation  in  a square  region  of  side  unity.  The  initial  data  for  the 
phase  field  and  composition  was  generated  from  a random  distribution  of  these  variables 
with  amplitude  0.01  about  a mean  value  of  a half.  The  following  values  were  used  for 
the  various  nondimensioncd  constants;  Pa  = — 0.1,^b  = O.I^Wa  = 10,  Wh  = 10,  Mi  = 
40,a7ide  = These  values  for  Mi  and  e are  chosen  to  He  in  the  asymptotic  regime  of  our 
analysis.  In  particular  the  size  of  e is  a compromise  between  making  it  sufficiently  small  that 
the  interfacial  layers  are  thin  enough  for  the  asymptotic  analysis  to  hold  but  large  enough 
that  they  are  adequately  resolved  by  the  finite  difference  grid.  In  Figure  8 we  display  the 
phase  field  at  times  t = 0,0.01,0.05,0.1  and  1. 

From  Figures  8(a-d),  it  is  clear  there  is  a rapid  initial  transient  evolution  of  the  phase 
field  to  ^ = 0 or  1,  which  is  complete  by  t = 0.1  when  fully  differentiated  regions  of  solid  and 
liquid  have  been  established.  This  is  in  agreement  with  our  asymptotic  analysis  although  it 
is  evident  that  there  is  also  a significant  coarsening  of  the  solid  regions  during  this  period. 
Subsequently  this  coctrsening  behaviour  continues  until  t = 1,  the  final  time  displayed. 
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We  monitored  the  total  amount  of  solute  present  as  a function  of  time.  It  was  found  that 
our  present  numerical  scheme  suffers  a small  loss  of  solute  from  the  system.  The  loss  is  most 
pronounced  at  small  times  when  the  system  is  evolving  most  rapidly.  A total  loss  of  1.26% 
by  t = 1 is  observed.  The  lack  of  solute  conservation  in  the  numerical  scheme  precludes  long- 
time integration  of  the  governing  equations,  since  the  loss  can  be  large  enough  to  eventually 
shift  the  average  solute  concentration  C to  values  outside  the  range  cs  < C < cl  for  the 
two-phase  system.  This  results  in  eventual  dissolution  of  one  of  the  phases,  leading  to  a 
final  state  consisting  of  a single  phase.  The  final  period  of  evolution  that  is  predicted  by  the 
asymptotic  analysis,  in  which  curvature  variations  are  important,  is  thus  difficult  to  observe 
with  the  current  numerical  scheme. 

The  qualitative  behaviour  observed  largely  substantiates  our  asymptotic  analysis,  demon- 
strating the  rapid  period  of  interfacial  genesis  and  order-one  motion  due  to  energy  differences. 
In  addition,  coarsening  is  observed  which  is  reminiscent  of  Ostwald  ripening. 

6 Conclusions 

In  this  paper  we  have  presented  a new  phase-field  model  which  describes  the  spatial  and 
temporal  evolution  of  isothermal  phase  transitions  between  ideal  binary  alloy  solution  phases. 
We  have  shown  how  the  parameters  which  characterise  this  phase-field  model  are  related  to 
material  parameters  and  thermodynamics  data  employed  in  classical  sharp  interface  models. 

We  conducted  an  asymptotic  analysis  of  the  governing  equations  in  the  distinguished 
limit  e — ► 0 with  Mi  = 0(e“^),  which  approximates  the  situation  when  the  mobility  of  the 
phase  field  is  large  and  its  diffusivity  small.  In  this  regime  there  exist  thin  interfacial  layers 
of  width  0{e)j  and  so  we  recovered  a sharp  interface  model  in  this  limit.  From  this  analysis 
we  determined  the  interfacial  dynamics  and  showed  that  the  interface  evolved  on  three  time 
scales.  There  is  short  initial  transient  in  which  interfaces  are  born.  It  results  primarily 
from  the  response  of  the  system  to  the  initial  data  we  have  chosen  and  does  not  necessarily 
represent  the  process  of  nucleation.  This  is  followed  by  evolution  on  the  diffusion  time  scale 
in  which  the  interfacial  motion  is  controlled  by  the  energy  difference  across  the  interface  and 
diffusion  in  the  solid  and  liquid  phases.  Finally  there  is  a long  period  in  which  the  interfaces 
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move  slowly  in  response  to  their  surface  energy  and  thermodynamic  disequihbrium.  This 
final  period  we  interpret  as  representing  Ostwald  ripening.  We  constructed  planar  interface 
solutions  in  this  asymptotic  hmit  and  were  able  to  characterise  the  interface  kinetics  in  terms 
of  the  thermodynamic  driving  force  and  show  a dependence  of  the  segregation  coefficient  on 
interfacial  velocity. 

Lastly,  we  presented  numerical  calculations  which  illustrate  the  predictions  of  our  asymp- 
totic analysis. 

It  is  clear  from  this  work  that  this  phase-field  model  provides  a good  description  of 
isothermal  phase  transitions  between  binary  solution  phases  and  provides  a first  step  in 
devising  a non-isothermal  phase-field  model  for  an  alloy  system. 
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8 Appendix:  The  Adjoint  Problem 

In  order  to  obtain  equation  (79),  representing  the  slow  time  evolution  of  a curved  interface, 
which  arises  from  solvabiHty  conditions  on  the  first  order  equations  (75)  and  (76)  we  require 
their  adjoint  solution  when  Vn  = 0.  Equations  (75)  and  (76)  when  K = 0 are: 

= , (97) 

+ = 0.  (98) 
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with 


^1)  ^ as  /9  ->  oo,  (99) 

^1)  1(^=0-)  as  ^ ^ — oo,  (100) 

We  introduce  the  adjoint  solution  [<j>'^ and  take  the  inner  product  of  (97)  and  (98)  with 
the  adjoint  solution.  This  gives 

J°°  Af-'  [^J°>  + 4>+dp  = 


+ 


f ~ + £'''[(^“^1  - ^°^)Cp)pfcc  - M*]}  dp+ 


(101) 


Thus  the  adjoint  problem  is: 


(c(«)(l  - = 0. 


This  has  a solution: 


p Bi  - B2c(°)(x) 

^0)(l  _ 

^+(/))  = B^mp), 


ip)  = / 

Jo 


dxj 


(102) 

(103) 

(104) 

(105) 


where  Bi  and  B2  are  constants  of  integration.  We  note  that  ~ a+p+^+H-o(l),  as  p ^ +00 
and  c'^  ~ a_p  + )0_  +0(1),  as  p — > —00.  The  constants  a+,  a_, /?+,/?_  are  determined  by  the 
constants  of  integration  as: 

Bi  B2 


a+  = 


a_  = 


C(°)  I (r=0+  ) ( 1 - c(°)  I (r=0+  ) ) 1 ” | (r=0+  ) 

Bi  ^2 

c(°)|(r=0-)(l  - c(°)|(r=o-))  1 “ c(°)  1(^=0") 

I3_  = Bi'y-  + B272", 


(106) 

(107) 

(108) 


13+  = Bi7i+  + ^272', 


(109) 
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where 


7i 


7r 


-£ 

-L 


c(0)(i  _ ^0))  c(0)  1(^=0-) (1  - c(°)|(r=o-)). 


^0)(1  - c(0))  c(0)|(r=o+)(l  - c(°)|(r=0+)). 


dp, 


dp, 


72 


72^ 


-L 

-I 


(l-c(0))  (1  - c(°)|(,=o-)). 


1 


(l-c(0))  (l-c(0)|(,=o+))J 

Inserting  these  forms  into  the  inner  product  (101)  gives 

cW(O-)  cW(0+) 


dp, 


dp. 


(110) 


(111) 


(112) 


(113) 


c(o)(o-)(l  - c(o)(0-))  cW(0+)(l  - c(°)(0+)) 

cW(0+)  c(i)(0-) 


(114) 


l_c(o)(0+)  l-c(“)(0-) 

= B2  r°°  M-i 

J 0 


(115) 


This  is  true  for  afl  choices  of  the  constants  Bi  and  B2  and  so  we  equate  the  coefficients  of 
these  quantities  in  the  above  expression  to  obtain  the  following  two  nontrivial  solvability 
conditions; 

c(i)(0-)  c(')(0+) 


c(o)(0-)(l  _ c(o)(0+))  c(o)(0+)(l  - c(o)(0-))’ 


(116) 


and 


:(')(0+) 


- 1 - }°)fo-i = /r 


1 - c(o)(0+)  1 - c(<>)(0-) 

These  solvability  conditions  may  be  combined  to  give: 


f-iSS— = -ro(^  + M-^Vn), 


(118) 


c(°)(0+)(l  - c(o)(0+)) 

where  Vn  = — t’t  is  the  the  slowly  varying  component  of  the  normal  velocity  of  the  interface, 
K,  = V^r  is  the  mean  curvature  of  the  interface,  and 


ro  = r 0°^fdp. 

J — 00 


(119) 
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Captions  to  the  Figures 


Figure  1.  The  free  energy  density  /(^,  T),  employed  by  Kobayashi  [16]  displayed  as  a func- 
tion of  (j)  for  three  values  of  the  temperature  near  the  melting  temperature,  Tm-  for  T < Tm 

{P  = -0.2),  T = Tm{P  = 0),  and  T > Tm  {P  = 0.2). 

Figure  2. The  free  energy  density  for  the  isothermal  solution  given  by  (10)  as  a function  of 
(f)  and  c. 

Figure  3.  (a)  The  common  tangent  construction  representing  equilibrium  of  a stationary 
flat  interface,  (b)  The  parallel  tangent  construction  predicted  by  this  model  for  a non  sta- 
tionary planar  interface.  In  both  cases  Wa  = Wb  = 100,  Pa  = —0.1,  and  Pb  = 0.1. 

Figure  4.  The  soHd  fraction,  So  versus  average  composition,  C,  obtained  from 

from  one- dimensional  stationary  solutions.  The  dashed  Hne  represents  two-phase  solutions. 
Figure  5.  The  solid  and  liquid  interfacial  concentrations  cs  and  cb  as  a function  of  growth 
rate,  V”,  for  a planar  constant  velocity  interface,  for  the  case  M = 1,  Wa  = 100,  Wb  = 100, 
Pa  — —0.1,  Pb  = 0.1.  The  dashed  curve  indicates  the  the  segregation  coefficient,  A;  as  a 
function  of  the  growth  rate. 

Figure  6.  The  one-dimensional  traveling  wave  solution  obtained  by  direct  numerical  inte- 
gration of  the  governing  equations.  The  concentration  profile  is  shown  for  t = 0, 1,  2,  3 (from 
left  to  right). 

Figure  7.  The  evolution  of  the  concentration  profile,  in  response  to  initial  conditions  in 
which  the  liquid  is  at  a spatial  uniform  concentration  away  from  its  equilibrium  value.  The 
concentration  profile  is  shown  dX  t — 0,2,4,6,8,10. 

Figure  8.  The  evolution  of  the  phase  field  (j)  in  time  and  space  on  the  domain  [0, 1]  x [0,1]. 
The  solid  curve  is  the  locus  of  points  ^ | and  represents  the  solid-liquid  interface.  The 

shaded  areas  represent  0<<^<^;^<(^<|;|<<^<|;|<(^<1  (in  increasing  order  of 
darkness).  Thus  the  darkest  regions  represent  solid  areas.  Plots  are  given  at  times  (a)  t = 0, 
(b)  t = 0.01,  (c)  t = 0.05,  (d)  t = 0.1,  (e)  t = 1. 
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Nickel 


(Jopp 


er 


Classical  Parameters 

Tm{K) 

L (J/cm^) 

Vjn  /mole)  “ 
a \ j / cm?)  ^ 
fi  (cm/Ks) 

D [cm? / s)  ® 

6 (cm) 


Phase  Field  Materials  Parameters 

e (J/cm)^ 

W (J/cm?) 

fi{T) 

Ml  (^cm^/Js) 

M2  (cm^/Js)  '* 


1728 

2350 

7.0 

3.7  X 10-5 
200  ' 

10-5 

5 X 10-®  f 


1358 

1725 

7.8 

2.8  X 10-5 
247'' 

10-5 

6.6  X 10-®  « 


4 X 10-5 

1.3  X 10^  7.2  X 105 

^ 'tB 


4.9  X 10« 
5.7  X 10-^ 


“ An  average  value  of  7.4  will  be  taken,  since  this  model  does  not  deal  with  volume  changes. 

^ Estimated  from  <7  = 0.7 L(vm/ Nq)^ , see  Coriell  and  Turnbull  [22];  Nq  is  Avagadro’s  Number. 
^ Estimated  from  fx  = ^^5^;  V,  is  the  speed  of  sound  in  liquid  metal 
(taken  as  3 x lO^cm/s);  Coriell  and  Turnbull  [22]. 

^ Calculated  from  using  equation  (40). 

* Typical  liquid  diffusion  coefficient. 

^ We  have  taken  the  interface  thickness  to  be  on  the  order  of  atomic  dimensions. 

^ Calculated  using  equation  (39). 

^ At  T*  = 154ZK. 


Table  1:  Material  and  phase-field  parameters  for  the  Ni-Cu  alloy. 
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